load packages

library(tidyverse)
library(knitr)
library(lme4)
library(lmerTest)

define color palettes and labels

pal_self_other = c("#FFA90A", "#247BA0")
pal_social_academic = c("#63647E", "#F25F5C")
pal_wave = c("#693668", "#A74482", "#F84AA7")
pal_label = c("#47A8BD", "#DBC057", "#FF3366")
pal_gender = c("#70c1b3","#247BA0")

parcel_labeller = labeller(label = c('social' = 'social parcels', 'other' = 'control parcels', 'self' = 'self parcels'),
                           domain = c('social' = 'social domain', 'academic' = 'academic domain'),
                           wave = c("t1" = "wave 1", "t2" = "wave 2", "t3" = "wave 3"))

label_df = expand.grid(label = c("social", "self", "other"),
              target = c("self", "other"),
              domain = c("social", "academic"),
              age = 13,
              expected_avg = 1,
              expected_diff = 1)

dcbw = #readRDS("~/dc_bw.Rds") +
  theme_classic() +
  theme(text = element_text(size = 14, family = "Futura Medium", color = "black"),
        panel.background = element_blank(),
        plot.background = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 14),
        legend.background = element_rect(fill = NA, color = NA),
        axis.line = element_line(color = "black"),
        axis.text = element_text(color = "black"),
        panel.grid.minor = element_blank())

load MRI data

  • exclude participants with >20% volumes with motion artifacts or low quality FX data
  • standardize betas within parcel (across participants and time)
  • recode standardized betas +/- 3 SDs from the mean as NA (0.8% of all observations)
# mri
# self_parcels = c(5, 17, 28, 47, 52, 54, 63, 66, 91, 116, 131, 139, 144, 147, 152, 156, 184, 198, 207, 225, 246, 249, 282, 292, 300, 301, 305, 309, 353, 354, 380, 399)
# social_parcels = c(4, 18, 20, 23, 29, 34, 49, 54, 59, 60, 62, 63, 66, 67, 76, 82, 111, 114, 116, 117, 139, 143, 146, 150, 152, 178, 179, 184, 189, 191, 203, 212, 224, 229, 236, 237, 238, 239, 241, 245, 246, 250, 259, 260, 266, 271, 273, 292, 301, 305, 309, 310, 322, 324, 328, 331, 333, 342, 343, 350, 353, 354, 369, 374, 391)

self_parcels = c(5, 17, 28, 47, 52, 66, 116, 147, 152, 156, 184, 198, 207, 225, 249, 292, 309, 354, 380)
social_parcels = c(18, 23, 29, 49, 54, 59, 62, 63, 67, 76, 111, 114, 139, 143, 146, 150, 178, 179, 189, 203, 212, 224, 229, 236, 238, 239, 245, 250, 259, 266, 271, 301, 305, 310, 322, 324, 328, 331, 333, 342, 343, 350, 374, 391)
                   
mri_exclusions = c('s002_t1', 's004_t1', 's008_t1', 's011_t1', 's017_t1', 
                   's026_t1', 's033_t2', 's034_t1', 's041_t1', 's044_t1', 
                   's047_t1', 's051_t1', 's054_t1', 's057_t1', 's059_t1', 
                   's061_t1', 's063_t1', 's070_t2', 's074_t1', 's074_t2', 
                   's078_t1', 's084_t1', 's090_t2', 's090_t3', 's094_t1', 
                   's094_t2', 's096_t1') 

parcellations = read_csv('../data/fxParcellations.csv') %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other'))) %>%
  mutate(wave = paste0("t", as.numeric(c(`10` = 1, `13` = 2, `16` = 3)[as.character(age)]))) %>%
  select(-age) %>%
  unite(sub_wave, c(subjectID, wave), remove = FALSE) %>%
  group_by(parcellation) %>%
  mutate(inclusion = ifelse(sub_wave %in% mri_exclusions, "excluded from MRI", "completed MRI"),
         beta = ifelse(sub_wave %in% mri_exclusions, NA, beta),
         sd = ifelse(sub_wave %in% mri_exclusions, NA, sd),
         beta_std = scale(beta, center = FALSE, scale = TRUE),
         mean_beta_std = mean(beta_std, na.rm = TRUE)) %>%
    select(-sub_wave) %>%
  ungroup()

# parcellations %>%
#   ggplot(aes(beta_std, fill = label)) +
#     geom_density(color = NA, alpha = .5) +
#     facet_grid(~wave) +
#     dcbw
# 
# psych::describe.by(parcellations$beta_std, group = parcellations$label)
# 
# parcellations %>%
#   filter(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3) %>%
#   group_by(subjectID, wave) %>%
#   summarize(n = n()) %>%
#   arrange(desc(n))
# 
# parcellations %>%
#   filter(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3) %>%
#   group_by(label) %>%
#   summarize(n = n()) %>%
#   ungroup() %>%
#   mutate(percent = round((sum(n) / nrow(filter(parcellations, !is.na(beta)))) * 100, 2))
# 
# parcellations %>%
#   filter(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3) %>%
#   group_by(label) %>%
#   summarize(meanOutlier = mean(beta_std, na.rm = TRUE),
#             sdOutlier = sd(beta_std, na.rm = TRUE))

parcellations_ex = parcellations %>%
  mutate(beta_std = ifelse(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3, NA, beta_std))

load demographic data and merge

# demographics
demo = read.csv("../data/SFIC_age.pds.gender.csv") %>%
  rename("subjectID" = SID,
         "wave" = wavenum,
         "gender" = Gender) %>%
  mutate(subjectID = sprintf("s%03d", subjectID),
         wave = paste0("t", wave),
         age_c = age - 13,
         age_c2 = age_c^2,
         pdss_c = pdss - 3,
         pdss_c2 = pdss_c^2)

# merge data
merged = parcellations_ex %>%
  full_join(., demo, by = c("subjectID", "wave")) %>%
  mutate(inclusion = ifelse(is.na(inclusion), "didn't complete MRI", inclusion)) %>%
  filter(!(subjectID == "s086" & wave == "t3")) #no MRI, task, or self-report data was collected

# subset data for modeling
neuro_model_data = merged %>%
  filter(!is.na(beta)) %>%
  select(subjectID, wave, age, age_c, age_c2, target, domain, parcellation, label, beta, beta_std)

neuro_model_data_ex = neuro_model_data  %>%
  na.omit()

# dummy code target and domain
neuro_model_data_dummy = neuro_model_data_ex %>%
  mutate(target = ifelse(target == "self", .5, -.5),
         domain = ifelse(domain == "social", .5, -.5))

sample descriptives

full sample

merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave) %>%
  summarize(n = n(),
            meanAge = mean(age, na.rm = TRUE),
            sdAge = sd(age, na.rm = TRUE)) %>%
  kable(format = 'pandoc', digits = 2)
wave n meanAge sdAge
t1 90 10.07 0.31
t2 57 13.04 0.31
t3 44 16.46 0.58
merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave, gender) %>%
  summarize(n = n()) %>%
  spread(gender, n) %>%
  kable(format = 'pandoc', digits = 2)
wave F M
t1 45 45
t2 30 27
t3 25 19

broken down by inclusion

age_table = merged %>% 
  select(subjectID, inclusion, wave, age, pdss, gender) %>%
  unique() %>%
  group_by(inclusion, wave) %>%
  summarize(N = n(),
            Mage = mean(age, na.rm = TRUE),
            SDage = sd(age, na.rm = TRUE),
            Mpdss = mean(pdss, na.rm = TRUE),
            SDpdss = sd(pdss, na.rm = TRUE))

gender_table = merged %>% 
  select(subjectID, inclusion, wave, age, pds, gender) %>%
  unique() %>%
  group_by(inclusion, wave, gender) %>%
  summarize(N = n()) %>%
  spread(gender, N) %>%
  rename("Female" = F,
         "Male" = M)

age_table %>%
  left_join(gender_table) %>%
  extract(wave, "Wave", "t([1-3]{1})") %>%
  arrange(Wave) %>%
  select(Wave, inclusion, N, Mage, SDage, Mpdss, SDpdss, Female, Male) %>%
  rename("MRI status" = inclusion) %>%
  kable(format = 'pandoc', digits = 2)
Wave MRI status N Mage SDage Mpdss SDpdss Female Male
1 completed MRI 57 10.08 0.32 2.16 0.74 33 24
1 didn’t complete MRI 12 10.00 0.25 2.00 0.77 4 8
1 excluded from MRI 21 10.07 0.33 1.83 0.58 8 13
2 completed MRI 44 13.06 0.33 3.57 1.05 25 19
2 didn’t complete MRI 8 13.07 0.20 3.69 0.65 4 4
2 excluded from MRI 5 12.79 0.28 2.62 0.85 1 4
3 completed MRI 34 16.33 0.46 4.75 0.39 19 15
3 didn’t complete MRI 9 17.06 0.60 4.78 0.67 6 3
3 excluded from MRI 1 15.72 NA 3.50 NA NA 1

plot

full sample

merged %>% 
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(19, 17, 21)) +
    labs(x = "\nage", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

MRI only

merged %>% 
  filter(!inclusion == "didn't complete MRI") %>%
  mutate(inclusion = ifelse(grepl("excluded", inclusion), "excluded", "included")) %>%
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(21, 19)) +
    labs(x = "\nage", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

descriptives

# collapsed across wave
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(domain, label) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
domain other M other SD self M self SD social M social SD
academic 0.01 0.92 -0.08 0.91 0.14 0.92
social -0.01 0.95 0.02 0.92 0.27 0.93
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(target, label) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
target other M other SD self M self SD social M social SD
other 0.00 0.94 -0.11 0.91 0.22 0.93
self -0.01 0.94 0.05 0.92 0.19 0.93
# as a function of wave
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(domain, label, wave) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
domain wave other M other SD self M self SD social M social SD
academic t1 0.00 0.92 -0.13 0.95 0.06 0.91
academic t2 0.00 0.96 -0.04 0.92 0.19 0.95
academic t3 0.03 0.89 -0.04 0.82 0.21 0.88
social t1 -0.02 0.94 0.00 0.91 0.22 0.91
social t2 0.00 0.96 0.00 0.95 0.30 0.98
social t3 -0.01 0.96 0.05 0.89 0.33 0.90
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(target, label, wave) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
target wave other M other SD self M self SD social M social SD
other t1 0.00 0.92 -0.10 0.92 0.16 0.90
other t2 0.00 0.97 -0.14 0.94 0.27 0.98
other t3 0.02 0.92 -0.10 0.85 0.27 0.88
self t1 -0.02 0.94 -0.03 0.94 0.11 0.92
self t2 0.00 0.94 0.11 0.92 0.22 0.95
self t3 0.01 0.93 0.11 0.86 0.27 0.90

visualize raw

hypothesis 0

domain_parc_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_raw = cowplot::plot_grid(domain_parc_plot_raw, target_parc_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, group = domain, color = domain)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = round(seq(-.4, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.45, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

soc_acad_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_social_academic[2], fill = pal_social_academic[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.4, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.45, .65)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

(h1_raw = cowplot::plot_grid(domain_plot_raw, soc_acad_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, group = target, color = target)) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.4, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

self_other_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(beta_std_avf = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, beta_std_avf) %>%
  unique() %>%
  spread(target, beta_std_avf) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_self_other[2], fill = pal_self_other[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.4, .5)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

(h2_raw = cowplot::plot_grid(target_plot_raw, self_other_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, beta_std, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = beta_std, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', linetype = '', fill = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, target, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_int_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, domain, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, target, beta_std_avg) %>%
  unique() %>%
  spread(target, beta_std_avg) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_raw = cowplot::plot_grid(int_plot_raw, soc_acad_int_plot_raw, self_other_int_plot_raw,
                       labels = c('A', 'B', 'C'), ncol = 3))

run models

domain x age

model_domain_equation = formula(beta_std ~ 1 + age_c*domain*label +
                                  age_c2*domain*label +
                                  (1 + age_c*domain || subjectID) +
                                  (1 + age_c*domain + age_c2*domain || parcellation))

model_domain_formula = lFormula(model_domain_equation, data = neuro_model_data_dummy)

model_domain_numFx = length(dimnames(model_domain_formula$X)[[2]])

model_domain_numRx = sum(as.numeric(lapply(model_domain_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_domain_maxfun = 10*(model_domain_numFx + model_domain_numRx + 1)^2

if (file.exists("../data/model_domain_dummy_2021.RDS")) {
  model_domain = readRDS("../data/model_domain_dummy_2021.RDS")
} else {
  model_domain = lmer(model_domain_equation, data = neuro_model_data_dummy, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_domain, "../data/model_domain_dummy_2021.RDS")
}

#summary(model_domain)

domain x target x age

model_target_equation = formula(beta_std ~ 1 + age_c*target*domain*label +
                                  age_c2*target*domain*label +
                                  (1 + age_c*target*domain || subjectID) +
                                  (1 + age_c*target*domain + age_c2*target*domain || parcellation))

model_target_formula = lFormula(model_target_equation, data = neuro_model_data_dummy)

model_target_numFx = length(dimnames(model_target_formula$X)[[2]])

model_target_numRx = sum(as.numeric(lapply(model_target_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_target_maxfun = 10*(model_target_numFx + model_target_numRx + 1)^2

if (file.exists("../data/model_target_dummy_2021.RDS")) {
  model_target = readRDS("../data/model_target_dummy_2021.RDS")
} else {
  model_target = lmer(model_target_equation, data = neuro_model_data_dummy, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_target, "../data/model_target_dummy_2021.RDS")
}
#summary(model_target)

compare models

anova(model_domain, model_target) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC X2 X2 df p
model_domain 29 462235.6 – – –
model_target 57 461963.7 327.95 28 < .001

summarize best fitting model

model_target %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept (age 13, label (control))", term),
         term = gsub("target", "Target", term),
         term = gsub("domain", "Domain", term),
         term = gsub("labelself", "Label (self)", term),
         term = gsub("labelsocial", "Label (social)", term),
         term = gsub("age_c", "Age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.3f [%.3f, %.3f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 3) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
Intercept (age 13, label (control)) -0.004 [-0.058, 0.049] 0.027 -0.157 389.202 .875
Age 0.003 [-0.003, 0.008] 0.003 0.975 165.622 .331
Target 0.003 [-0.014, 0.021] 0.009 0.364 179.471 .716
Domain 0.009 [-0.016, 0.033] 0.013 0.690 144.735 .491
Label (self) -0.034 [-0.242, 0.175] 0.107 -0.315 356.532 .753
Label (social) 0.242 [0.100, 0.385] 0.073 3.328 356.480 .001
Age2 0.000 [-0.001, 0.001] 0.001 0.280 599.799 .780
Age x Target 0.002 [-0.003, 0.007] 0.002 0.748 73.548 .457
Age x Domain -0.002 [-0.009, 0.005] 0.004 -0.503 63.368 .617
Target x Domain -0.009 [-0.044, 0.026] 0.018 -0.527 164.410 .599
Age x Label (self) 0.007 [-0.008, 0.022] 0.008 0.924 353.133 .356
Age x Label (social) 0.018 [0.008, 0.029] 0.005 3.511 352.651 .001
Target x Label (self) 0.255 [0.198, 0.313] 0.029 8.664 2045.761 < .001
Target x Label (social) -0.051 [-0.090, -0.012] 0.020 -2.540 2030.155 .011
Domain x Label (self) 0.066 [-0.003, 0.135] 0.035 1.877 1009.675 .061
Domain x Label (social) 0.118 [0.071, 0.165] 0.024 4.908 1004.290 < .001
Target x Age2 -0.002 [-0.004, -0.000] 0.001 -2.108 1427.460 .035
Domain x Age2 -0.003 [-0.005, -0.001] 0.001 -3.046 1884.460 .002
Label (self) x Age2 0.001 [-0.003, 0.006] 0.002 0.600 358.197 .549
Label (social) x Age2 -0.004 [-0.007, -0.001] 0.001 -2.925 356.212 .004
Age x Target x Domain 0.001 [-0.009, 0.012] 0.005 0.263 69.501 .794
Age x Target x Label (self) 0.026 [0.013, 0.039] 0.007 3.911 186305.317 < .001
Age x Target x Label (social) 0.006 [-0.003, 0.015] 0.005 1.393 186299.729 .164
Age x Domain x Label (self) -0.006 [-0.019, 0.007] 0.007 -0.914 186304.467 .361
Age x Domain x Label (social) -0.005 [-0.013, 0.004] 0.005 -1.016 186297.730 .310
Target x Domain x Label (self) -0.022 [-0.132, 0.088] 0.056 -0.393 186311.567 .694
Target x Domain x Label (social) 0.097 [0.022, 0.172] 0.038 2.547 186303.865 .011
Target x Domain x Age2 0.000 [-0.004, 0.004] 0.002 0.057 1492.001 .955
Target x Label (self) x Age2 -0.011 [-0.017, -0.004] 0.003 -3.164 186307.919 .002
Target x Label (social) x Age2 0.004 [-0.000, 0.009] 0.002 1.852 186304.113 .064
Domain x Label (self) x Age2 0.006 [-0.000, 0.013] 0.003 1.909 186314.601 .056
Domain x Label (social) x Age2 0.005 [0.001, 0.010] 0.002 2.372 186302.120 .018
Age x Target x Domain x Label (self) -0.006 [-0.032, 0.020] 0.013 -0.421 186304.080 .673
Age x Target x Domain x Label (social) -0.014 [-0.031, 0.004] 0.009 -1.507 186298.997 .132
Target x Domain x Label (self) x Age2 0.009 [-0.005, 0.022] 0.007 1.257 186307.980 .209
Target x Domain x Label (social) x Age2 -0.001 [-0.010, 0.008] 0.005 -0.288 186301.661 .773
# summarize random effects
print(VarCorr(model_target), comp = c("Variance","Std.Dev."), digits = 3)
##  Groups          Name                 Variance Std.Dev.
##  parcellation    target:domain:age_c2 0.000000 0.00000 
##  parcellation.1  age_c:target:domain  0.000000 0.00000 
##  parcellation.2  domain:age_c2        0.000000 0.00000 
##  parcellation.3  target:age_c2        0.000000 0.00000 
##  parcellation.4  target:domain        0.000000 0.00000 
##  parcellation.5  age_c:domain         0.000000 0.00000 
##  parcellation.6  age_c:target         0.000000 0.00000 
##  parcellation.7  age_c2               0.000032 0.00566 
##  parcellation.8  domain               0.008080 0.08989 
##  parcellation.9  target               0.001437 0.03791 
##  parcellation.10 age_c                0.000859 0.02931 
##  parcellation.11 (Intercept)          0.198861 0.44594 
##  subjectID       age_c:target:domain  0.000829 0.02878 
##  subjectID.1     target:domain        0.005907 0.07686 
##  subjectID.2     age_c:domain         0.000553 0.02351 
##  subjectID.3     age_c:target         0.000139 0.01181 
##  subjectID.4     domain               0.004338 0.06586 
##  subjectID.5     target               0.001406 0.03749 
##  subjectID.6     age_c                0.000158 0.01256 
##  subjectID.7     (Intercept)          0.002164 0.04651 
##  Residual                             0.665033 0.81550

simple slopes

run models

# self social > academic
self_social = emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c2", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "quadratic"))

# social self > other
social_self = emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c2", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "quadratic"))

make table

social_self %>%
  bind_rows(self_social) %>%
  select(contrast, parcel, age_effect, estimate, SE, df, t.ratio, p.value) %>%
  rename("b" = estimate,
         "t" = t.ratio,
         "p" = p.value) %>%
  mutate(b = round(b, 3) * -1, #flip signs for it's .5 - (-.5)
         SE = round(SE, 3),
         df = round(df, 2),
         t = abs(round(t, 2)),
         p = round(p, 3)) %>%
  kable(format = "pandoc")
contrast parcel age_effect b SE df t p
social self > other self linear 0.026 0.010 6704.64 2.70 0.007
social self > other self quadratic -0.008 0.005 133230.55 1.79 0.073
self social > academic social linear -0.013 0.007 684.47 1.76 0.079
self social > academic social quadratic 0.002 0.003 57999.28 0.53 0.595

visualize fitted models

reForm = as.formula("~(1 + age_c*target*domain + age_c2*target*domain || parcellation)")

neuro_plot_data = with(neuro_model_data_dummy, 
                    expand.grid(target = unique(target), 
                                domain = unique(domain),
                                parcellation = unique(parcellation),
                                age = unique(age),
                                stringsAsFactors = F)) %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other')),
         age_c = age - 13, 
         age_c2 = age_c^2, 
         subjectID = NA)
neuro_plot_data$expected = predict(model_target, newdata = neuro_plot_data, re.form = reForm)
neuro_plot_data$expected_mean = predict(model_target, newdata = neuro_plot_data, re.form = NA)

neuro_plot_data = neuro_plot_data %>%
  mutate(target = factor(target, levels = c(-.5, .5), labels = c("other", "self")),
         domain = factor(domain, levels = c(-.5, .5), labels = c("academic", "social")))

hypothesis 0

domain_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(name = "", values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n", color = "") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_fitted = cowplot::plot_grid(domain_parc_plot, target_parc_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label, color = NA), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .4)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_social_academic[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .4)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h1_fitted = cowplot::plot_grid(domain_plot, soc_acad_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, target, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label, color = NA), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .3)) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_self_other[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .3)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h2_fitted = cowplot::plot_grid(target_plot, self_other_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_mean, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = 1, se = FALSE, size = 1.25) + 
  #geom_point() + 
  scale_y_continuous(breaks = c(-.4, .2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '', linetype = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, .2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))


self_other_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, .2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_fitted = cowplot::plot_grid(int_plot, soc_acad_int_plot, self_other_int_plot,
                       labels = c('A', 'B', 'C'), ncol = 3))

compare raw and fitted models

hypothesis 0

raw

h0_raw

fitted

h0_fitted

hypothesis 1

raw

h1_raw

fitted

h1_fitted

hypothesis 2

raw

h2_raw

fitted

h2_fitted

hypothesis 3

raw

h3_raw

fitted

h3_fitted